#include "./rand.h"
#include <math.h>

/**
 http://school.anhb.uwa.edu.au/personalpages/kwessen/shared/Marsaglia99.html
 */

#if defined(__INTEL_COMPILER)
#	pragma warning ( disable : 981 )
#endif

uint32_t mds_hash32( uint32_t a)
{
	a = (a+0x7ed55d16) + (a<<12);
	a = (a^0xc761c23c) ^ (a>>19);
	a = (a+0x165667b1) + (a<<5);
	a = (a+0xd3a2646c) ^ (a<<9);
	a = (a+0xfd7046c5) + (a<<3);
	a = (a^0xb55a4f09) ^ (a>>16);
	return a;
}

typedef
struct
{
	uint32_t z;
	uint32_t w;
	uint32_t jsr;
	uint32_t jcong;
} mds_rand_state;

static mds_rand_state mds_rand_ = 
{
	/* z= */     362436069,
	/* w= */     521288629, 
	/* jsr= */   123456789, 
	/* jcong= */ 380116160
};


static
uint32_t mds_MWC( mds_rand_state *s )
{
	s->z = 36969 * (s->z & 65535) + (s->z >> 16);
    s->w = 18000 * (s->w & 65535) + (s->w >> 16);
    return (s->z << 16) + s->w;  /* 32-bit result */
}

static uint32_t mds_SHR3( mds_rand_state *s )
{
	s->jsr^=(s->jsr<<17);
	s->jsr^=(s->jsr>>13);
	s->jsr^=(s->jsr<< 5);
	return s->jsr;
}


static uint32_t mds_CONG( mds_rand_state *s)
{
	s->jcong=69069* s->jcong+1234567;
	return s->jcong;
}

static uint32_t mds_KISS( mds_rand_state *s )
{
	return ( mds_MWC(s) ^ mds_CONG(s) ) + mds_SHR3(s);
}


static void mds_INIT( mds_rand_state *s, uint32_t seed )
{
	s->w     = ( seed = mds_hash32(seed) );
	s->z     = ( seed = mds_hash32(seed) );
	s->jsr   = ( seed = mds_hash32(seed) );
	s->jcong =  mds_hash32(seed);
}

#if 0
size_t mds_rand_state_length(void)
{
	return sizeof( mds_rand_state );
}

void  *mds_rand_state_buffer(void)
{
	return & mds_rand_;
}
#endif

void mds_srand( uint32_t seed )
{
	mds_INIT( &mds_rand_ , seed );
	
}

uint32_t mds_rand(void)
{
	return mds_KISS( &mds_rand_ );
}

Real mds_alea(void)
{
	static const Real fac= 1.0 / 4294967296.0;
	return ( (Real)0.5 + (Real)mds_rand() ) * fac;
}

Real mds_normal(void)
{
	static const Real two_pi   = 6.28318530717959;
	static int        generate = 1;
	static Real       g2       = 0;
	if( generate != 0)
	{
		const Real x1     = mds_alea();
		const Real x2     = mds_alea();
		const Real angle  = two_pi * x2;
		const Real half_r = -(Real)log( x1 );
		const Real radius = (Real)sqrt( half_r + half_r );
		const Real g1     = radius * (Real)cos( angle );
		g2                = radius * (Real)sin( angle );
		generate = 0;
		return g1;
	}
	else 
	{
		generate = 1;
		return g2;
	}
	
}
